dpath<-"/home1/30/jc227089/LD_inv/ran_outs"

setwd("/home1/30/jc227089/LD_inv/results")


# for plotting mean LD at N and D loci on the invasion front against neigh.size and generation
#takes dframe out of LD_sims.r
front.wframe<-function(dframe, percent=F){
	require(lattice)
	D.data<-tapply(dframe$D.out, list(dframe$gen, dframe$neigh.size), mean, na.rm=T)
	N.data<-tapply(dframe$N.out, list(dframe$gen, dframe$neigh.size), mean, na.rm=T)
	diff.data<-D.data-N.data
	if (percent==T) diff.data<-100*diff.data/N.data
	intl<-length(N.data[,1])%/%length(N.data[1,])
	x<-seq(from=1, length.out=length(N.data[1,]), by=intl)
	x<-matrix(x, ncol=length(x), nrow=length(N.data[,1]), byrow=T)
	y<-matrix(1:length(N.data[,1]), ncol=length(x[1,]), nrow=length(N.data[,1]))
	z<-diff.data
	filled.contour(x=seq(from=1, length.out=length(N.data[1,]), by=intl), y=1:length(diff.data[,1]), t(z), 
		zlim=c(0, max(z)), color.palette=gray.colors)
}

# Takes dframe out of LD_sims.R
#plots histograms of D' for DD vs NN
front.hist<-function(dframe){
	dc.lev<-levels(as.factor(dframe$dc))
	gen.lev<-levels(as.factor(dframe$gen))
	for (dd in 1:length(dc.lev)){
		for (gg in 1:length(gen.lev)){
				temp<-subset(dframe, dframe$dc==dc.lev[dd] & dframe$gen==gen.lev[gg])
				if (length(temp$R)==0) next
				xmax<-max(temp$R, na.rm=T)
				Xbreaks<-seq(0, xmax, length.out=20)
				ymax<-max(c(hist(temp$R[which(temp$p.type=="NN")])$density, hist(temp$R[which(temp$p.type!="NN")])$density), na.rm=T)
				fig.id<-paste("DhistDC=", 10*as.numeric(dc.lev[dd]), "GEN=", gen.lev[gg], ".png", sep="")
				png(fig.id, width=15, height=15, units="cm", res=300)
				par(cex.lab=1.5, mar=c(5,5,2,2))
				hist(temp$R[which(temp$p.type=="NN")], col="grey30", breaks=Xbreaks, xlab="D'", main=fig.id, ylim=c(0,ymax), freq=F)
				par(new=T)
				over.col<-col2rgb("grey70") #convert colour to rgb with 50% transparency
				over.col<-rgb(red=over.col[1], green=over.col[1], blue=over.col[1], alpha=125, maxColorValue=255)
				hist(temp$R[which(temp$p.type!="NN")], breaks=Xbreaks, col=over.col, ylab="", xlab="", main="", axes=F, ylim=c(0,ymax), freq=F)
				legend(0.5*xmax, 0.9*ymax, bty="n", legend=c("NN", "DD"), title="Locus pair type", fill=c("grey30", over.col))
				dev.off()		
		}	
	}	
}

#takes ibd dframe from LD_sims and plots IBD slopes for each locus
ibd.comp<-function(dframe){
	require(reshape)
	dframe<-melt(dframe, id.vars=1:5)
	dc.lev<-levels(as.factor(dframe$dc))
	gen.lev<-levels(as.factor(dframe$gen))
	for (dd in 1:length(dc.lev)){
		for (gg in 1:length(gen.lev)){
			temp<-subset(dframe, dframe$dc==dc.lev[dd] & dframe$gen==gen.lev[gg])
			if (length(temp$variable)==0) next
			fig.id<-paste("IbyDDC=", 10*as.numeric(dc.lev[dd]), "GEN=", gen.lev[gg], ".png", sep="")
			png(fig.id, width=15, height=15, units="cm", res=300)
			par(cex.lab=1.5, mar=c(5,5,2,2))
			plot(temp$variable, temp$value)
			dev.off()
		}
	}		
}

# takes ld dframe and summarizes LD for each locus across each rep/gen/dc
ld.comp<-function(dframe, nloci){
	dc.lev<-levels(as.factor(dframe$dc))
	gen.lev<-levels(as.factor(dframe$gen))
	nreps<-max(dframe$rep)
	dframe$rep<-as.factor(dframe$rep)
	out<-c()
	for (dd in 1:length(dc.lev)){
		for (gg in 1:length(gen.lev)){
			for (ll in 1:nloci){
				temp<-subset(dframe, dframe$dc==dc.lev[dd] & dframe$gen==gen.lev[gg] & (dframe$loc1==ll | dframe$loc2==ll))
				byrep<-tapply(temp$R, temp$rep, mean)
				byrep<-data.frame(dc=dc.lev[dd], gen=gen.lev[gg], rep=1:nreps, locus=ll, mean.R=byrep)
				out<-rbind(out, byrep)
			}
		}
	}
	return(out)	
}



fles<-list.files(dpath)
ibd<-fles[grep("IbyD", fles)]
ldf<-fles[grep("ld", fles)]

conc<-c()
for (ii in 1:length(ldf)){
	load(paste(dpath, "/ldout", ii, ".RData", sep=""))
	ld.out$rep<-ld.out$rep+3*(ii-1)
	conc<-rbind(conc, ld.out)
	rm(ld.out)	
}
ld.reps<-ld.comp(conc, 20)
ld.reps$dc<-as.numeric(as.character(ld.reps$dc))
ld.reps$gen<-as.numeric(as.character(ld.reps$gen))
ld.reps$locus<-as.numeric(as.character(ld.reps$locus))

conc<-c()
for (ii in 1:length(ibd)){
	load(paste(dpath, "/IbyDout", ii, ".RData", sep=""))
	IbyD.out$rep<-IbyD.out$rep+3*(ii-1)
	conc<-rbind(conc, IbyD.out)
	rm(IbyD.out)	
}
ibd.reps<-melt(conc, id.vars=1:5)
names(ibd.reps)[6:7]<-c("locus", "mean.b")
ibd.reps$locus<-as.numeric(ibd.reps$locus)

joined<-merge(ld.reps, ibd.reps)
save(joined, file="LDplusIBD.RData")



joined.100<-subset(joined, joined$gen==100)
plot(joined.100$mean.R, joined.100$mean.b)
points(joined.100$mean.R[which(joined.100$locus>10)], joined.100$mean.b[which(joined.100$locus>10)], col="red")

for (ii in 1:30){
	rws<-which(joined.100$rep==ii)
	temp<-joined.100[rws,]
	png(paste("plot", ii, ".png", sep=""), height=15, width=15, units="cm", res=300)
	plot(temp$mean.R, temp$mean.b)
	points(temp$mean.R[which(temp$locus>10)], temp$mean.b[which(temp$locus>10)], col="red")
	dev.off()	
}